PRO PLOTS_BGPS_GALPOS_COMP
  
  ;;==========================================================================
  ;; Routine to make comparison plots for galactic positions derived
  ;;   from the NEAR/FAR determination of the code.  This analysis
  ;;   routine is a heavily-modified version of INVESTIGATE_GALPOS.pro
  
  ;; Define the COMMON BLOCK
  COMMON FFORE_BLOCK,n,rb3,R0,d,R,Z,tau,f_data,sig_f,corr,farlist,do_Tfit,$
     lpstr,rho_hi,rho_h2,rho_star
  COMMON DATA_BLOCK, s,nh2_irdc,nh2_bgps,sig_nb,sig_ni
  
  ;; Since we are using routines buried in BGPS_FFORE, compile now to
  ;; they are available.
  FORWARD_FUNCTION BGPS_FFORE_MODEL, BGPS_FFORE_LOOPIE
  RESOLVE_ROUTINE, 'BGPS_FFORE', /COMPILE_FULL_FILE, /EITHER
  
  ;; Here's where you can change the SUBSET used (i.e. TEMP,
  ;; training set, etc.)  ALSO the SUFF is which rotation curve used
  subset = ''
  
  do_Tfit = 0b
  
  IF n_elements(rb3) EQ 0 THEN $
     restore,'irdc_dist_model/bgps_rb3'+subset+'.sav',/ver
  n = n_elements(rb3)
  
  ;; Define arrays for TEMP_DEPENDENT_VARIABLES
  tau    = dblarr(n)
  f_data = dblarr(n)
  sig_f  = dblarr(n)
  
  ;; Set up some of the variables needed throughout the program
  chi    = cgSymbol('chi',/PS)
  gr_rho = cgSymbol('rho',/PS)
  gr_tau = cgSymbol('tau',/PS)
  
  ;; Get galactic params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Using !MW parameters -- Goes into COMMON block!
  R0 = !MW.R0
  d = dindgen(!MW.NBINS*5.)*!MW.BINSIZE + !MW.BINSTART
  corr = R0 / 8.5d3
  
  ;;==========================================================================
  ;; Set up the plotting environment   
  myps,'./irdc_dist_model/analysis_plots/bgps_galpos_plot_comp'+subset+'.eps',$
       ct=0,xsize=7.5,ysize=6,/cmyk
  
  multiplot,[2,1],mxtitle='Galactocentric Distance [kpc]',$
            mytitle='Galactocentric Distance [kpc]'
  
  FOR rotc=0b, 3b, 2b DO BEGIN
     
     ;; Do re-usable calculations
     GENERATE_LPSTR, ROTC=rotc, SUFF=suff, RTIT=rtit
     
     ;; There are currently 3 cases of the model
     nmods     = 3
     dfn       = strarr(nmods)
     thisnf    = bytarr(n)
     fitdist   = fltarr(n)
     ffore_vec = fltarr(n)
     
     cnames = ['YGB7','Lime Green','TG3','RED7']
     colors = cgColor(cnames)
     textcolor = cgColor('RYB8')
     thick=[6,5,4,3]
     
     tstring = 'FIT'
     
     message,'Dust temp = '+tstring,/inf
     
     ytit=''                    ;'Galactocentric Distance [kpc]'
     ytf = (rotc EQ 0) ? "(I0)" : ""
     
     plotsym,0,2.0                    ;,/fill
     plot,[-9.0,1.0],[-9,7],/nodata,$ ;xtit='Galactocentric Distance [kpc]',$
          ytit=ytit,ytickformat=ytf,/xst,/yst,xticks=10,xticklen=0.005,$
          tit=rtit,charsize=1.2,/isotropic,xtickformat='funky'
  
     
     chmin = fltarr(nmods)
     ndof  = fltarr(nmods)
     tdp   = fltarr(nmods)
     
     ;; Solar Circle and Tangent Circles
     ells = findgen(1001)/1.e3*180.
     scl = findgen(1001)/1.e3*360.
     xells = -(R0/1.d3 * cos(ells*!dtor)) *  sin(ells * !dtor)
     yells = + (R0/1.d3 * cos(ells*!dtor)) * cos(ells * !dtor) - R0/1.d3
     xsc = -R0/1.d3 * cos(scl * !dtor)
     ysc = -R0/1.d3 * sin(scl * !dtor)
     oplot,xells,yells,linestyle=1,color=cgColor('BLK6'),thick=1
     oplot,xsc,ysc,linestyle=1,color=cgColor('BLK6'),thick=1
     
     ;; l=30 line (through end of the Bar)
     plots,0,-R0/1.d3
     plots,-R0/1.d3*cos(30*!dtor),R0/1.d3*sin(30*!dtor),/cont,linestyle=2,$
           thick=3,color=cgColor('BLK5')
     
     ;; FOR j=6, nmods-1 DO BEGIN
     j = 1                      ; Model Case 2 is where it's at!
     
     message,'Case #'+string(j+1,format="(I0)"),/inf
     
     dfn[j] = './irdc_dist_model/data/case'+$
              string(j+1,format="(I0)")+$
              '_irdc_TdFIT'+subset+suff+'.sav'
     
     restore,dfn[j]
     chmin[j] = bestnorm
     ndof[j] = dof
     
     message,'Minimum Chisq = '+string(bestnorm,format="(F0.2)"),/inf
     
     do_Tfit = 1b
     
     residj = BGPS_FFORE_LOOPIE( res, NEARFAR=thisnf , FFORE_VEC=ffore_vec)
     
     ;; Sort out near vs. far vectors, and derive fitdist
     far = thisnf & near = ~ thisnf
     fitdist = near * lpstr.kdn + far * lpstr.kdf
     
     ;; Plot only for well-constrained sources!
     cind    = where(constrain,nc)
     rb3     = rb3[cind]
     fitdist = fitdist[cind]
     print,'Plotting for '+string(nc,format="(I0)")+' well-constrained sources'
     
     ;; Make x,y position for plotting
     x = -fitdist/1.d3 * sin(rb3.l*!dtor)
     y = -R0/1.d3 + fitdist/1.d3 * cos(rb3.l*!dtor)
     FOR kk = 0L, nc-1 DO $
        plots,x[kk],y[kk],color=cgColor('Black'),psym=8,symsize=rb3[kk].c_meas
     
     
     ;; Sun and Galactic Center
     plots,0,0,psym=7,color=cgColor('BLK5'),symsize=1.2,thick=8
     plots,0,-R0/1.d3,psym=7,color=cgColor('BLK5'),symsize=1.2,thick=8
     cgText,0.1,0.4,'GC',charsize=0.7
     cgText,0.1,-8.1,'Sun',charsize=0.7
     
     cgText,-2.3,2.0,'Sagittarius Arm?',charsize=0.7,$
            alignment=0.5,orientation=40.,color='Red'
     cgText,-5.0,3.3,'Perseus Arm?',charsize=0.7,$
            alignment=0.5,orientation=30.,color='Red'
     cgText,-7.5,5.5,'Outer Arm?',charsize=0.7,$
            alignment=0.5,orientation=30.,color='Red'
     cgText,-1.0,-3.0,'End of Bar?',charsize=0.7,$
            alignment=0.5,orientation=0.,color='Red'
     
     plotsym,0,0.5,/fill
     
     multiplot
     
  ENDFOR
  
  
  myps,/done,/mp
  
END
